#include <iostream>
#include <math.h>
#include <sys/stat.h>

#include <vtkDICOMImageReader.h>
#include <vtkImageData.h>
#include <vtkSmartPointer.h>
#include <vtkMarchingContourFilter.h>
#include <vtkDecimatePro.h>
#include <vtkPolyData.h>
#include <vtkPolyDataNormals.h>
#include <vtkPolyDataMapper.h>
#include <vtkSTLWriter.h>
#include <vtkAppendPolyData.h>
#include <vtkSmoothPolyDataFilter.h>
#include <vtkQuadricClustering.h>

using namespace std;

int otsu(vtkImageData * volumeData,bool bkgIgnored=false);

int otsu(vtkImageData * volumeData,bool bkgIgnored){
	cout << "otsu start------" << endl;
	int *dims=volumeData->GetDimensions();
    int width=dims[0];
    int height=dims[1];
    int frame=dims[2];
    cout << "width:" << width << ",height:" << height << ",frame:" << frame <<endl;
    int componentSum=width*height*frame;
    cout << "componentSum:" << componentSum <<endl;
    double  *scalarRange=volumeData->GetScalarRange();
    cout << "scalarRange=[" << scalarRange[0] <<","<<scalarRange[1]<<"]"<<endl;
    int level = static_cast<int>(scalarRange[1]-scalarRange[0])+1;
    cout << "level:" << level <<endl;
    int threshold = static_cast<int>(scalarRange[0]);
    float w0,w1,u0,u1,u0tmp,u1tmp,deltaTmp,deltaMax=0;
    int *componentCount=new int[componentSum];
    memset(componentCount,0,sizeof(int)*componentSum);
    float *componentPercent=new float [componentSum];
    memset(componentPercent,0.0,sizeof(float)*componentSum);
    short *volumeDataPtr= static_cast<short *>(volumeData->GetScalarPointer());
    cout << "volumeDataPtr:" << volumeDataPtr <<endl;
    int diff=-static_cast<int>(scalarRange[0]);
    cout << "diff:" << diff <<endl;
    for (int i = 0; i <frame ; ++i) {
        for (int j = 0; j <height ; ++j) {
            for (int k = 0; k <width ; ++k) {
                int index=i*width*height+j*width+k;
                short value=volumeDataPtr[index];
                componentCount[value+diff]++;
            }
        }
    }
    if(bkgIgnored){
        componentCount[0]=0;
    }
    for (int i = 0; i < level; ++i) {
        componentPercent[i]=componentCount[i]*1.0f/componentSum;
    }
    for (int i = 0; i <level ; ++i) {
        w0=w1=u0=u1=u0tmp=u1tmp=deltaTmp=0.0;
        for (int j = 0; j <level; ++j) {
            //Foreground
            if(j<=i){
                w0+=componentPercent[j];
                u0tmp+=j*componentPercent[j];
            }
            //Background
            else{
                w1+=componentPercent[j];
                u1tmp+=j*componentPercent[j];
            }
            u0=u0tmp/w0;
            u1=u1tmp/w1;
            deltaTmp= static_cast<float>(w0*w1*pow((u0-u1),2));
            if(deltaTmp>deltaMax){
                deltaMax=deltaTmp;
                threshold= static_cast<int>(i-diff);
            }
        }
	}
	delete [] componentCount;
    delete [] componentPercent;
	return threshold;
}

int dcm2stl(const char *dcmPath,const char *stlName)
{
    
    cout << "Start convert dicoms to stl from：" << dcmPath <<endl;
    cout << "Check file is exits" <<endl;
    if (FILE *file = fopen(dcmPath, "r")) {
        cout << "File exits:"  << dcmPath << endl;
        fclose(file);
    } else {
        cout << "No found file:"  << dcmPath << endl;
        return -1;
    } 
    try{
        vtkSmartPointer<vtkDICOMImageReader> dcmReader = vtkSmartPointer<vtkDICOMImageReader>::New();
        dcmReader->SetDataByteOrderToLittleEndian();
        dcmReader->SetDirectoryName(dcmPath);
        dcmReader->SetDataSpacing(3.2,3.2,1.5);
        dcmReader->Update();
    	cout << "dcmReader->Update()------" << endl;
    	
         vtkSmartPointer<vtkImageData> volumeData = dcmReader->GetOutput();
		int threshold=otsu(volumeData);
		double  *scalarRange=volumeData->GetScalarRange();
		int *dims=volumeData->GetDimensions();
		int dataSize=dims[0]*dims[1]*dims[2];
		short *volumeDataPtr= static_cast<short *>(volumeData->GetScalarPointer());
		for(int i=0;i<dataSize;++i){
			if(volumeDataPtr[i]<threshold){
				volumeDataPtr[i]=scalarRange[0];
			}
		}
		threshold=otsu(volumeData,true);
		cout << "threshold=" << threshold <<endl;
		vtkSmartPointer<vtkMarchingContourFilter> surfaceExtractor = vtkSmartPointer<vtkMarchingContourFilter>::New();
		surfaceExtractor->SetInputData(volumeData);
		surfaceExtractor->SetValue(0,threshold);
		surfaceExtractor->Update();
		cout << "surfaceExtractor->Update()------" << endl;

		vtkSmartPointer<vtkSmoothPolyDataFilter> surfaceSmoother = vtkSmartPointer<vtkSmoothPolyDataFilter>::New();
		surfaceSmoother->SetInputData(surfaceExtractor->GetOutput());
		surfaceSmoother->SetNumberOfIterations(100);
		surfaceSmoother->Update();
		cout << "surfaceSmoother->Update()------" << endl;

		vtkSmartPointer<vtkPolyDataNormals> surfaceNormals = vtkSmartPointer<vtkPolyDataNormals>::New();
		surfaceNormals->SetInputData(surfaceSmoother->GetOutput());
		surfaceNormals->FlipNormalsOn();
		surfaceNormals->Update();

		vtkSmartPointer<vtkQuadricClustering> deciFilter = vtkSmartPointer<vtkQuadricClustering>::New();
		deciFilter->SetInputData(surfaceNormals->GetOutput());
		deciFilter->SetNumberOfXDivisions(200);
		deciFilter->SetNumberOfYDivisions(200);
		deciFilter->SetNumberOfZDivisions(200);
		deciFilter->Update();
		vtkSmartPointer<vtkPolyDataMapper> boneMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		boneMapper->SetInputData(deciFilter->GetOutput());
		boneMapper->ScalarVisibilityOff();
		boneMapper->Update();
		cout << "boneMapper->Update()------" << endl;
		vtkSmartPointer<vtkSTLWriter>  stlWriter = vtkSmartPointer<vtkSTLWriter>::New();
		vtkSmartPointer<vtkAppendPolyData> apd= vtkSmartPointer<vtkAppendPolyData>::New();
		apd->AddInputData(boneMapper->GetInput());
		apd->Update();
		cout << "vtkAppendPolyData->Update()------" << endl;
		stlWriter->SetInputData(apd->GetOutput());
		stlWriter->SetFileName(stlName);
		stlWriter->SetFileTypeToBinary();
		stlWriter->Write();
		cout << "Success to write  stl: " << stlName <<endl;
		}catch(...){
			cout << "Failed to convert stl"<<endl;
			return -1;
		}													
    return 0;
}






